//      nmf_bcd.cc
//      Block Coordinate Descent NMF algorithm
//      
//      Copyright 2012 tqlong <tqlong@espada>
//      
//      This program is free software; you can redistribute it and/or modify
//      it under the terms of the GNU General Public License as published by
//      the Free Software Foundation; either version 2 of the License, or
//      (at your option) any later version.
//      
//      This program is distributed in the hope that it will be useful,
//      but WITHOUT ANY WARRANTY; without even the implied warranty of
//      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//      GNU General Public License for more details.
//      
//      You should have received a copy of the GNU General Public License
//      along with this program; if not, write to the Free Software
//      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
//      MA 02110-1301, USA.

#include "helper/stl_helper.h"
#include "helper/arma_helper.h"
#include "helper/mpi_helper.h"

#include "nmf.h"

extern po::variables_map vm;
extern mpi::communicator* p_world;

void nmf_bcd()
{
  RequiredParameter(vm, "bcd")("input")("output")("log")("k")("alpha")("iter");
  reset_seed();

  std::string inFile = VM(vm, "input", std::string);
  std::string logFile = VM(vm, "log", std::string);
  matrix_t A;
  read_sparse_matrix<matrix_t>(inFile, A, vm.count("binary"));

  int m = A.n_rows, n = A.n_cols, k = VM(vm, "k", int);
  double alpha = VM(vm, "alpha", double);
  int maxiter = VM(vm, "iter", int);

  matrix_t W, H; 
  gen_sparse_matrix<matrix_t>(W, m, k, 0.5);
  gen_sparse_matrix<matrix_t>(H, k, n, 0.5);

  matrix_t A_H(A.n_rows, H.n_rows);
  matrix_t A_W(A.n_cols, W.n_cols);
  matrix_t W_T_W(W.n_cols, W.n_cols);
  matrix_t H_T_H(H.n_rows, H.n_rows);
  matrix_t Wnew = W, Hnew = H;
  
  FileAccessCreate(logFile);
  for (int iter = 0; iter < maxiter; iter++) {
    W = nmf_bcd_W<matrix_t>(A, W, H, alpha, A_H, H_T_H, Wnew);
    H = nmf_bcd_H<matrix_t>(A, W, H, alpha, A_W, W_T_W, Hnew);
    
    double err = arma::norm(A-W*H,"fro");
    FileAccessAppend(log, logFile);
    log << "\t" << iter << "\t" << err << "\t" << TOTAL_SECONDS << ENDL;
    if (iter % 10 == 0 || iter == maxiter-1)
      MSG("iter=" << iter << " err=" << err << " d=" << duration_all());
  }
}
